/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 Mario Orsi
   This file is part of Brahms.
   Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
   Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
   You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. */

/************************************************************************************************************
 * lpp.c -- Direct evaluation of the intramembrane pressure ("lateral pressure profile") 
 ************************************************************************************************************/

#include "dataStructs.h"

extern const VecR *siteCoords; 
extern const VecR *siteVelocity; 
extern const Site *site; 
extern const VecR region; 
extern const int nSites, stepCount, sizeHistLpp, limitLpp, nDOPCsDSPCs, nDOPEs, contourIrvKirk,  
  halfSizeHistLpp; // defined in startUp.c
extern const real deltaT, *mass, timeNow, lppDeltaZ, extTemperature;
extern int countLpp; // global, defined & initialised in startUp.c
extern real *totalHistLpp, *histLppTau1, *histLppTau2, *histPress_xx, *histPress_yy, *histPress_zz, 
  meanLppSlabVol, meanLppDeltaZ, meanLppBoxHeight; // global, defined and initialised in startUp.c

void ZeroLppHistograms(); // external function declaration (function defined in startUp.c)

/**************************************************
 * DistributeVirial -- local stress decomposition *
 **************************************************/
void DistributeVirial( const RMat virTens, const int i, const int j ) 
{       
  int iSlab, jSlab;
  real virial_xx = virTens.u[0]; 
  real virial_yy = virTens.u[4]; 
  real virial_zz = virTens.u[8]; 
  real latPress = 0.5*(virial_xx + virial_yy) - virial_zz; // pL - pN

  iSlab = siteCoords[ i ].z / lppDeltaZ + halfSizeHistLpp; // make all slab bins positive
  jSlab = siteCoords[ j ].z / lppDeltaZ + halfSizeHistLpp; 

  if (contourIrvKirk) { // use the Irving-Kirkwood contour

    int rightmostSlab, leftmostSlab, nInterSlabs, slabCounter;
    real zij, fLatPress;
    real weightFac; // weighting factor for distributing the virial
    VecR rij;
    VSub( rij, siteCoords[ i ], siteCoords[ j ] );
    VWrap( rij, z );   /* Minimum Image Convention */
    zij = fabs( rij.z );

    if ( iSlab == jSlab ) {/* if both sites are in the same slab */
      totalHistLpp[ iSlab ] += latPress; /* place the entire virial in the common slab */
      histPress_xx[ iSlab ] += virial_xx;
      histPress_yy[ iSlab ] += virial_yy;
      histPress_zz[ iSlab ] += virial_zz;
    } else {
      rightmostSlab = Max( iSlab, jSlab );
      leftmostSlab = Min( iSlab, jSlab ); 
      nInterSlabs = ( rightmostSlab - leftmostSlab ) - 1; // compute the #slabs between the sites, without considering MIC/PBC
      if ( nInterSlabs == 0 || nInterSlabs == (sizeHistLpp-2) ) { // if directly adjacent, or through PBC
	weightFac = 0.5; // assign half virial to each of the two slabs
	totalHistLpp[ iSlab ] += weightFac * latPress; 
	totalHistLpp[ jSlab ] += weightFac * latPress;
	histPress_xx[ iSlab ] += weightFac * virial_xx;
	histPress_xx[ jSlab ] += weightFac * virial_xx;
	histPress_yy[ iSlab ] += weightFac * virial_yy;
	histPress_yy[ jSlab ] += weightFac * virial_yy;
	histPress_zz[ iSlab ] += weightFac * virial_zz;
	histPress_zz[ jSlab ] += weightFac * virial_zz;
      } else { // this is the general case when the two sites are separated by one or more slabs
	if ( nInterSlabs < halfSizeHistLpp ) { // if PBC are definitely not involved  
	  weightFac = 1.0 / (nInterSlabs+2);  
	  totalHistLpp[ iSlab ] += weightFac * latPress; // in primis, assign the adequate bit of virial to iSlab and jSlab
	  totalHistLpp[ jSlab ] += weightFac * latPress; 
	  histPress_xx[ iSlab ] += weightFac * virial_xx;
	  histPress_xx[ jSlab ] += weightFac * virial_xx;
	  histPress_yy[ iSlab ] += weightFac * virial_yy;
	  histPress_yy[ jSlab ] += weightFac * virial_yy;
	  histPress_zz[ iSlab ] += weightFac * virial_zz;
	  histPress_zz[ jSlab ] += weightFac * virial_zz;
	  for ( slabCounter = leftmostSlab + 1; slabCounter < rightmostSlab; slabCounter++ ) { // loop through the inter-slabs between i and j
	    totalHistLpp[ slabCounter ] += weightFac * latPress;
	    histPress_xx[ slabCounter ] += weightFac * virial_xx;
	    histPress_yy[ slabCounter ] += weightFac * virial_yy;
	    histPress_zz[ slabCounter ] += weightFac * virial_zz;
	  }
	} else { // PBC wrapping 
	  weightFac = 1.0 / (sizeHistLpp-nInterSlabs);
	  totalHistLpp[ iSlab ] += weightFac * latPress; 
	  totalHistLpp[ jSlab ] += weightFac * latPress; 
	  histPress_xx[ iSlab ] += weightFac * virial_xx;
	  histPress_xx[ jSlab ] += weightFac * virial_xx;
	  histPress_yy[ iSlab ] += weightFac * virial_yy;
	  histPress_yy[ jSlab ] += weightFac * virial_yy;
	  histPress_zz[ iSlab ] += weightFac * virial_zz;
	  histPress_zz[ jSlab ] += weightFac * virial_zz;
	  for ( slabCounter = 0; slabCounter < leftmostSlab; slabCounter++ ) { // scan the bottom slabs 
	    totalHistLpp[ slabCounter ] += weightFac * latPress;
	    histPress_xx[ slabCounter ] += weightFac * virial_xx;
	    histPress_yy[ slabCounter ] += weightFac * virial_yy;
	    histPress_zz[ slabCounter ] += weightFac * virial_zz;
	  }
	  for ( slabCounter = rightmostSlab + 1; slabCounter < sizeHistLpp; slabCounter++ ) { // scan the top slabs
	    totalHistLpp[ slabCounter ] += weightFac * latPress;
	    histPress_xx[ slabCounter ] += weightFac * virial_xx;
	    histPress_yy[ slabCounter ] += weightFac * virial_yy;
	    histPress_zz[ slabCounter ] += weightFac * virial_zz;
	  }
	}
      } 
    } 
  } else { /* use the Harasima contour */
    totalHistLpp[ iSlab ] += .5 * latPress;
    totalHistLpp[ jSlab ] += .5 * latPress;
    histPress_xx[ iSlab ] += .5 * virial_xx;
    histPress_xx[ jSlab ] += .5 * virial_xx;
    histPress_yy[ iSlab ] += .5 * virial_yy;
    histPress_yy[ jSlab ] += .5 * virial_yy;
    histPress_zz[ iSlab ] += .5 * virial_zz;
    histPress_zz[ jSlab ] += .5 * virial_zz;
  }
}

static void PrintLpp()
{
  int bin;
  real distanceFromBilayerCentre;
  FILE *lppFilePtr;
  char lppFileName[ 20 ];
  char time_ns[ 9 ]; // string containing timeNow, e.g. "100", in ns

  // Print lateral pressure profile lpp = PL - PN = 0.5*(Pxx+Pyy) - Pzz 
  sprintf( time_ns, "%d", ( int ) ( 0.5 + timeNow *  TIME_ns ) );
  lppFileName[ 0 ] = '\0';
  strcat( lppFileName, "lpp-" );
  strcat( lppFileName, time_ns );
  strcat( lppFileName, ".dat" );
  lppFilePtr = fopen( lppFileName, "w" );
  
  for ( bin = 0; bin < sizeHistLpp; bin++ ) {
    distanceFromBilayerCentre = ( bin + 0.5 ) * meanLppDeltaZ - meanLppBoxHeight / 2.; 
    fprintf( lppFilePtr, "%7.3f", ( ( distanceFromBilayerCentre ) * LENGTH_nm ) );
    fprintf( lppFilePtr, "%15.3f", ( totalHistLpp[ bin ] * PRESSURE_atm ) ); 
    fprintf( lppFilePtr, "\n" ); 
  }
  fclose( lppFilePtr );

  // Print Pxx(z)
  sprintf( time_ns, "%d", ( int ) ( 0.5 + timeNow *  TIME_ns ) );
  lppFileName[ 0 ] = '\0';
  strcat( lppFileName, "Pxx-" );
  strcat( lppFileName, time_ns );
  strcat( lppFileName, ".dat" );
  lppFilePtr = fopen( lppFileName, "w" );
  
  for ( bin = 0; bin < sizeHistLpp; bin++ ) {
    distanceFromBilayerCentre = ( bin + 0.5 ) * meanLppDeltaZ - meanLppBoxHeight / 2.;
    fprintf( lppFilePtr, "%7.3f", ( ( distanceFromBilayerCentre ) * LENGTH_nm ) );
    fprintf( lppFilePtr, "%15.3f", ( histPress_xx[ bin ] * PRESSURE_atm ) );
    fprintf( lppFilePtr, "\n" );
  }
  fclose( lppFilePtr );

  // Print Pyy(z)
  sprintf( time_ns, "%d", ( int ) ( 0.5 + timeNow *  TIME_ns ) );
  lppFileName[ 0 ] = '\0';
  strcat( lppFileName, "Pyy-" );
  strcat( lppFileName, time_ns );
  strcat( lppFileName, ".dat" );
  lppFilePtr = fopen( lppFileName, "w" );
  
  for ( bin = 0; bin < sizeHistLpp; bin++ ) {
    distanceFromBilayerCentre = ( bin + 0.5 ) * meanLppDeltaZ - meanLppBoxHeight / 2.;
    fprintf( lppFilePtr, "%7.3f", ( ( distanceFromBilayerCentre ) * LENGTH_nm ) );
    fprintf( lppFilePtr, "%15.3f", ( histPress_yy[ bin ] * PRESSURE_atm ) );
    fprintf( lppFilePtr, "\n" );
  }
  fclose( lppFilePtr );

  // Print Pzz(z)
  sprintf( time_ns, "%d", ( int ) ( 0.5 + timeNow *  TIME_ns ) );
  lppFileName[ 0 ] = '\0';
  strcat( lppFileName, "Pzz-" );
  strcat( lppFileName, time_ns );
  strcat( lppFileName, ".dat" );
  lppFilePtr = fopen( lppFileName, "w" );
  
  for ( bin = 0; bin < sizeHistLpp; bin++ ) {
    distanceFromBilayerCentre = ( bin + 0.5 ) * meanLppDeltaZ - meanLppBoxHeight / 2.;
    fprintf( lppFilePtr, "%7.3f", ( ( distanceFromBilayerCentre ) * LENGTH_nm ) );
    fprintf( lppFilePtr, "%15.3f", ( histPress_zz[ bin ] * PRESSURE_atm ) );
    fprintf( lppFilePtr, "\n" );
  }
  fclose( lppFilePtr );

  sprintf( time_ns, "%d", ( int ) ( 0.5 + timeNow *  TIME_ns ) );
  lppFileName[ 0 ] = '\0';
  strcat( lppFileName, "Tau1-" );
  strcat( lppFileName, time_ns );
  strcat( lppFileName, ".dat" );
  lppFilePtr = fopen( lppFileName, "w" );
  
  for ( bin = 0; bin < sizeHistLpp; bin++ ) {
    distanceFromBilayerCentre = ( bin + 0.5 ) * meanLppDeltaZ - meanLppBoxHeight / 2.; 
    fprintf( lppFilePtr, "%7.3f", ( ( distanceFromBilayerCentre ) * LENGTH_nm ) );
    fprintf( lppFilePtr, "%15.3f", ( histLppTau1[ bin ] * PRESSURE_atm * LENGTH_nm ) ); 
    fprintf( lppFilePtr, "\n" ); 
  }
  fclose( lppFilePtr );

  sprintf( time_ns, "%d", ( int ) ( 0.5 + timeNow *  TIME_ns ) );
  lppFileName[ 0 ] = '\0';
  strcat( lppFileName, "Tau2-" );
  strcat( lppFileName, time_ns );
  strcat( lppFileName, ".dat" );
  lppFilePtr = fopen( lppFileName, "w" );
  
  for ( bin = 0; bin < sizeHistLpp; bin++ ) {
    distanceFromBilayerCentre = ( bin + 0.5 ) * meanLppDeltaZ - meanLppBoxHeight / 2.; 
    fprintf( lppFilePtr, "%7.3f", ( ( distanceFromBilayerCentre ) * LENGTH_nm ) );
    fprintf( lppFilePtr, "%15.3f", ( histLppTau2[ bin ] * PRESSURE_atm * Sqr( LENGTH_nm ) ) ); 
    fprintf( lppFilePtr, "\n" ); 
  }
  fclose( lppFilePtr );
}

static void ComputeSurfaceTension() { // integrate the lateral pressure profile: gamma = - 0.5 * int( lpp ) dz
  int bin;
  FILE *gammaFilePtr;
  real gammaLeft = 0.;  // surface tension of "left" monolayer
  real gammaRight = 0.; // surface tension of "right" monolayer
  
  for ( bin = halfSizeHistLpp - 1; bin >= 0; bin-- )  // "left" monolayer - going backwards from bilayer centre to water
    gammaLeft  -= totalHistLpp[ bin ];   // accumulate pN-pL = -(pL-pN) = -pi
  
  for ( bin = halfSizeHistLpp; bin < sizeHistLpp; bin++ ) // "right" monolayer
    gammaRight -= totalHistLpp[ bin ];   // accumulate pN-pL = -(pL-pN) = -pi
    
  gammaLeft *= meanLppDeltaZ; // normalise
  gammaRight *= meanLppDeltaZ; // normalise

  gammaFilePtr = fopen( "surfaceTensionPerMonolayer.log", "a");
  fprintf( gammaFilePtr, "Surface tension of \"left\"  monolayer from integrating lpp: gammaLeft  = %10.6f dyn/cm\n",   ( gammaLeft  * TENSION_dyn_cm ) );
  fprintf( gammaFilePtr, "Surface tension of \"right\" monolayer from integrating lpp: gammaRight = %10.6f dyn/cm\n\n", ( gammaRight * TENSION_dyn_cm ) );
  fclose( gammaFilePtr );
} 

/********************************************************************************************************************************************
 * ComputeLatPressMoments -- computes the first & second moment of the lateral pressure profile [Cantor, Biophys J, 80, 2284 (2001)]
 ********************************************************************************************************************************************/
static void ComputeLatPressMoments() 
{
  const real T_room_K = 293;  // room temperature [K], same as 20 Celsius
  const real kBT_room_J = kB_IN_J_OVER_K * T_room_K; // thermal energy at room temperature ~ 4.0453E-21 J

  const real DOPC_BILAYER_BENDING_MODULUS_J = 0.76e-19;  // [J] from Pan et al., BJ [2008]
  const real DOPE_BILAYER_BENDING_MODULUS_J = 0.90e-19;  // [J] from Chen & Rand, BJ [1997]
  /*   const real c0_DOPC_per_nm = -0.115; // [1/nm] DOPC spontaneous curvature, from Chen & Rand, BJ, [1997] */
  /*   const real tau1m_DOPC_kBT_per_nm = 0.5*DOPC_BILAYER_BENDING_MODULUS_J/kBT_room_J*c0_DOPC_per_nm; // pressure profile's first moment */

  real p, tau1m, tau2m, kG_left, kG_right, kG, tau1mLeft, tau1mRight, tau2mLeft, tau2mRight, z;
  real c0_left, c0_right, c0; // monolayer spontaneous curvatures: left monolayer, right monolayer, average

  real xi = 1.; // [nm] z-distance of the pivotal plane (polar-apolar interface) from the bilayer centre [Orsi et al, JPCM, 2010]
  real kBendingMonolayer_J; /* monolayer bending rigidity modulus in units of J */

  int bin;
  FILE *fPtr;
  
  tau1mLeft = tau1mRight = tau2mLeft = tau2mRight = 0.;

  // Setting the *monolayer* bending modulus for subsequent c0 calculation - this formula gives a rough estimate for mixtures too
  kBendingMonolayer_J = 0.5 * ( nDOPCsDSPCs * DOPC_BILAYER_BENDING_MODULUS_J + nDOPEs * DOPE_BILAYER_BENDING_MODULUS_J ) / ( nDOPCsDSPCs + nDOPEs );
  //  printf ("***Just checking - kBendingMonolayer_J = %e\n",kBendingMonolayer_J );
  
  for ( bin = halfSizeHistLpp - 1; bin >= 0; bin-- ) { // "left" monolayer - going backwards from bilayer centre to water
    z = ( ( halfSizeHistLpp - 1 ) - bin + 0.5 ) * meanLppDeltaZ;  // compute the (positive!) distance "z" from the bilayer centre
    //    printf("\nz = %f\n", z * LENGTH_nm);
    p = totalHistLpp[ bin ];   // get the local net lateral pressure acting at distance "z"
    histLppTau1[ bin ] = - z * p; 
    tau1mLeft -= histLppTau1[ bin ]; //accumulate first integral moment of the lateral pressure (negative summation as we are on the negative side of the z-axis)
    histLppTau2[ bin ] = Sqr( z ) * p;
    tau2mLeft += histLppTau2[ bin ]; // accumulate second integral moment of the lateral pressure
  }
  
  tau1mLeft *= meanLppDeltaZ;
  tau2mLeft *= meanLppDeltaZ;

  c0_left = tau1mLeft / ( kBendingMonolayer_J / ENERGY_J ) / LENGTH_nm; // units: [ 1 / nm ] 
  kG_left = 2. * xi * tau1mLeft - tau2mLeft;
  
  fPtr = fopen( "lppMoments.log", "a");
  /*   fprintf( fPtr, "tau1m left layer = % 7.4f kB * T / nm\t", tau1mLeft * ENERGY_J / kBT_room_J ); */
  /*   fprintf( fPtr, "tau2m left layer = % 7.4f kB * T\n",     tau2mLeft * ENERGY_J / kBT_room_J ); */
  /*   fprintf( fPtr, "c0_left = tau1m_left / k = % 5.4f / nm\t", c0_left ); */
  /*   fprintf( fPtr, "kG_left = 2*xi*tau1m - tau2m_left = % 7.4f kBT\n\n", kG_left * ENERGY_J / kBT_room_J ); */

  for ( bin = halfSizeHistLpp; bin < sizeHistLpp; bin++ ) { // "right" monolayer
    z = ( bin - halfSizeHistLpp + 0.5 ) * meanLppDeltaZ;  // compute the distance from the bilayer centre
    p = totalHistLpp[ bin ];   // get the local net lateral pressure acting at distance "z"
    histLppTau1[ bin ] = z * p;
    tau1mRight += histLppTau1[ bin ];
    histLppTau2[ bin ] = Sqr( z ) * p;
    tau2mRight += histLppTau2[ bin ];
  }

  tau1mRight *= meanLppDeltaZ;
  tau2mRight *= meanLppDeltaZ;
 
  c0_right = tau1mRight / ( kBendingMonolayer_J / ENERGY_J ) / LENGTH_nm; // units: [ 1 / nm ] 
  kG_right = 2. * xi * tau1mRight - tau2mRight;

  /*   fprintf( fPtr, "tau1m right layer = % 7.4f kB * T / nm\t", tau1mRight * ENERGY_J / kBT_room_J ); */
  /*   fprintf( fPtr, "tau2m right layer = % 7.4f kB * T\n",     tau2mRight * ENERGY_J / kBT_room_J ); */
  /*   fprintf( fPtr, "c0_right = tau1m_right / k = % 5.4f / nm\t", c0_right ); */
  /*   fprintf( fPtr, "kG_right = - ( tau2m_right - 2*xi*tau1m_right ) = % 7.4f kBT\n\n", kG_right * ENERGY_J / kBT_room_J ); */

  tau1m = 0.5 * ( tau1mRight + tau1mLeft );
  tau2m = 0.5 * ( tau2mRight + tau2mLeft );
  c0 = tau1m / ( kBendingMonolayer_J / ENERGY_J ); // units: [ 1 / nm ] 
  kG = 2. * xi * tau1m - tau2m;

  fprintf( fPtr, "Monolayer torque tension tau = tau1m = % .2f kB * T_room / nm", tau1m * ENERGY_J / kBT_room_J );
  if ( nDOPCsDSPCs ) fprintf( fPtr, " - for DOPC this should range from -0.09 to -2.39 kBT/nm" );
  fprintf( fPtr, "\nMonolayer spontaneous curvature c0 = tau1m / k = % .2f / nm (assuming k = %.2e J)\n", c0, kBendingMonolayer_J );
  fprintf( fPtr, "Monolayer Gaussian curvature kG=2*xi*tau1m-tau2m=%.2f kBT (-%.2f kBT < kG < 0) ", kG*ENERGY_J/kBT_room_J, kBendingMonolayer_J/kBT_room_J ); 
  fprintf( fPtr, "(assuming xi = %.2f nm)\n", xi ); 
  fclose( fPtr );

  fPtr = fopen( "tau1m.dat", "a");
  fprintf( fPtr, "% 8.6f\n", tau1m * ENERGY_J / kBT_room_J );  // kB * T_room / nm
  fclose( fPtr );

  fPtr = fopen( "tau2m.dat", "a");
  fprintf( fPtr, "% 8.6f\n", tau2m * ENERGY_J / kBT_room_J ); //  kB * T
  fclose( fPtr );

} 

void EvalLpp() 
{
  RMat kinEnDensTensor; // kinetic energy tensor density tensor - used to in the evaluation (accumulation) of the pressure profile kinetic contribution 
  int n, bin; // "bin" is the histogram bin number

  real slabVolume = region.x * region.y * lppDeltaZ;   // computing the volume over which the density is worked out
  
  DO_SITE { // computing and sorting kin en contribution to (pL - pN)
    VOuter( kinEnDensTensor.u, siteVelocity[ n ], siteVelocity[ n ] );
    MScale( kinEnDensTensor.u, mass[ site[ n ].type ] );
    bin = siteCoords[ n ].z / lppDeltaZ + halfSizeHistLpp;    
    totalHistLpp[ bin ] += 0.5 * (kinEnDensTensor.u[0]+kinEnDensTensor.u[4]) - kinEnDensTensor.u[8];
    histPress_xx[ bin ] += kinEnDensTensor.u[0];
    histPress_yy[ bin ] += kinEnDensTensor.u[4];
    histPress_zz[ bin ] += kinEnDensTensor.u[8];
  }
  
  meanLppSlabVol += slabVolume; // accumulating - really only necessary for NPT (redundant at constant volume, but its no big deal 2 compute, hence leave it)
  meanLppDeltaZ += lppDeltaZ;  // lppDeltaZ computed in SingleStep()
  meanLppBoxHeight += region.z; 
  
  ++countLpp;

  if ( limitLpp == countLpp ) {  // if happy with the amount of data collected
    for ( bin = 0; bin < sizeHistLpp; bin++ )  {// normalising to have an averaged value
      totalHistLpp[ bin ] /= meanLppSlabVol; 
      histPress_xx[ bin ]  /= meanLppSlabVol; 
      histPress_yy[ bin ]  /= meanLppSlabVol; 
      histPress_zz[ bin ]  /= meanLppSlabVol; 
    }
    
    meanLppDeltaZ /= countLpp;
    meanLppBoxHeight /= countLpp; 
    
    ComputeSurfaceTension();
    ComputeLatPressMoments();

    PrintLpp();

    countLpp = 0;
    meanLppSlabVol = 0.;
    ZeroLppHistograms();
  }
}
